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I. Introduction 


Recent advances in supersonic and hypersonic aerospace technology have led to a renewed 
interest in the stability and transition to turbulence of high speed flows. The last 30 years 
have intermittedly witnessed some vigorous attempts to understand some of the funda- 
mental routes to transition for incompressible flows. While a fairly comprehensive picture 
of the initial stages leading to the breakdown of an incompressible laminar boundary layer 
has emerged (mostly under controlled conditions) [9], the non-linear effects responsible for 
transition at high speeds are still very much a mystery. However, the current nonlinear 
incompressible theories, numerical simulations and experiments will, hopefully, serve as a 
guide in gaining a better understanding of the mechanisms present in the supersonic and 
hypersonic regimes. 

Lees and Lin [13] were the first to study the linear stability of inviscid compressible 
flows in the 1940’s. Their study was restricted to low Mach numbers and the conclusions 
drawn only apply to the first unstable mode. Mack unraveled some of the intricacies among 
the many possible modes present when the full stability equations are considered (includ- 
ing viscous and conductivity terms). From 1960 to the present time, Mack has devoted 
his energies to a thorough numerical investigation of the linear stability characteristics of 
inviscid compressible flows at both subsonic and supersonic Mach numbers. An extensive 
review of compressible linear theory appears in his 1984 AGARD report [16]. His research 
led to the discovery of the higher unstable modes for Mach numbers above 2.2. The impor- 
tance of the higher modes stems from their high spatial and temporal amplification rates 
which make them susceptible to noise from the free-stream or from windtunnel walls. One 
of the main differences between first and second mode waves is their response to cooling. 
Wall cooling stabilizes first mode perturbation waves, while second mode perturbations 
are destabilized. At Mach numbers greater than 7, cone experiments have confirmed that 
the second mode is primarily responsible for the observed instabilities [22,3]. Therefore, 
new techniques must be developed to control the growth rates in these regimes. 

The existence of an inflectional point in the mean boundary layer profile, a necessary 
and sufficient condition to insure the instability of inviscid incompressible flows, remains 
valid at all Mach numbers if the relative Mach number is subsonic everywhere *. For 
compressible flows, the inflection point is replaced by the generalized inflection point which 
takes into account the mean density variation across the boundary-layer. However, there 
are also non-inflectional waves which depend on the existence of local pockets of supersonic 
relative Mach number in the boundary-layer. These waves can be unstable and can exist 
at all supersonic free-stream Mach numbers. Multiple modes also exist for these non- 
inflectional waves. The above linear stability characteristics of compressible flows are fully 
described in Mack’s review articles [16,17]. 

I'The relative Mach number measures the difference between the mean velocity in the direction of the 
wave vector and the phase velocity with respect to the sound speed. 
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To date, there are only a handful of experiments that have studied the detailed linear 
stability characteristics of flat plates [12] and cones [10]. Laufer and Vrebalovich studied the 
linear stability of a Mach 2.2 flow on a flat plate, and theirs is perhaps the most complete 
(if not the only) flat plate experiment to date at supersonic velocities. Not only are hot 
wire measurements difficult to perform and to interpret correctly [19], but the stability 
characteristics of the boundary layer are affected by free-stream disturbances emanating 
from the wind tunnel walls. Since the disturbances scale as M it becomes increasingly 
difficult to obtain good, quantitative data at the higher speeds. The quiet Mach 6 tunnel, 
currently under construction at NASA Langley Research Center [l] should alleviate some 
of these concerns. 

Among the few numerical codes geared towards the analysis of compressible stability, 
the best known axe perhaps those of Mack [15] and Malik [18]. The former calculates 
spatial and temporal amplification rates, while the latter is strictly restricted to temporal 
analysis; a Gaster transformation is then applied to calculate the temporal and spatial 
growth rates. A new spectral stability code has also been developed at NASA Langley [14]. 
Until the work by Erlebacher and Hussaini [6,7], there have been no non-linear simulations. 
Erlebacher and Hussaini studied the evolution of two 3-D waves superimposed on a 2-D 
wave. These were first mode waves. They numerically demonstrated the existence of a 
secondary instability at Mach 4.5. The current incompressible weakly non-linear theory of 
Herbert [8] has been extended to include compressible flow by El-Hady [4] and Nayfeh [20]. 
Nayfeh treats both subsonic and supersonic regimes, while El-Hady only presents results 
up to Mach 1.2. Both studies are based on a quasi-parallel theory and treat subharmonic 
disturbances. 

In this paper, a recent code developed to solve numerically the non-stationary, three- 
dimensional compressible Navier-Stokes equations applied to wall-bounded flows, is used 
to study the nonlinear evolution of one two-dimensional second mode at Mach 4.5. A 
possible application of perturbation analysis to this problem is also discussed. 


II. Equations 


The compressible Navier-Stokes equations for the stability of parallel flows are, in dimen- 
sionless form, (Erlebacher and Hussaini [6]) 

+ V-(pv) = 0 (1) 

+ V.(pvv) = (2) 

% + VP + 7PV - v = — V-(mVT) + (7 - 1)$ (3) 


where 


t = p(Vv + Vv T ) — -^(V-v) I 


(4) 
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is the viscous stress and 


( 5 ) 


$ = I(VtT+ VtT r ) : 7 

is the viscous dissipation function. The equation of state is 

iMlp = pT, (6) 

with 7 = C p /C v , the ratio of specific heats. At high Mach number, temperature variations 
across the boundary layer become important and the temperature dependence of viscosity 
must be taken into account. The Prandtl number is assumed equal to 0.7, and Suther- 
land’s law is prescribed for viscosity. Lengths are non-dimensionalized with respect to the 
displacement thickness 6*. Velocity, temperature, viscosity, and density are normalized 
with respect to free-stream values. Pressure is normalized with respect to the total head 

In incompressible boundary-layer flows, it has been observed that the evolution to 
laminar breakdown for ribbon-excited Tollmien-Schlichting (TS) waves occupies a space 
of several TS wavelengths [11]. Apart from numerical considerations, a direct simulation 
of a 3-D boundary layer flow far into the non-linear regime is not currently practical, 
even on today’s most powerful computers. The situation is worse for compressible flows 
where the instability mechanisms are attenuated. One possible approximation is to neglect 
the growth of the mean boundary-layer. The perturbation flow is then periodic and the 
expense of the simulation is greatly decreased. This has been successful in incompressible 
theories [8] and numerical simulations [24]. One of the major reasons for this success is that 
the actual transition process occurs within a very short streamwise extent. At high Mach 
numbers, the growth rates are attenuated, the transition region is longer, and it is not clear 
whether this approximation is valid. Nonetheless, in an initial attempt to gain insight into 
the physics of second mode wave evolution, a temporal simulation is considered. The 
hope is that the observed features are at least qualitatively correct. The computational 
domain is restricted to one period in the streamwise distance. Periodicity is also assumed 
in the spanwise direction. A properly posed initial-value problem is set up, and qualitative 
comparisons are made with experimental data by relating time (in the simulations) to the 
streamwise direction (x) in the experiments. In order that the mean parallel boundary 
layer flow, (p m (y), U m (y),T m (y)), be a stationary solution to the Navier-Stokes equations, 
source terms must be added to equations (2)-(3). The x momentum equation must be 
supplemented by the term 



while the forcing term 


must be added to the energy equation (3) . In obtaining these relations, the mean pressure 
is kept constant across the boundary layer. The mean velocity has only a streamwise 
component. 



7 = 1 r. 

p RePrM dy 
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III. Algorithm 


The algorithm described below, has been implemented in 3-D although only 2-D results 
are presented in this paper. Three-dimensional simulations of first mode wave interactions 
have been detailed elsewhere [7]. The system of equations (1) -(3) is solved in conservative 
form under the parallel flow assumption. Periodic boundary conditions in the streamwise 
and spanwise directions permit a Fourier representation of the primitive variables v, p and 
p. For example, 

N z /2 N ,/ 2 

v(x,y,z,t) = £ £ (9) 

k=-N z / 2+1 l=-N z / 2+1 

where N x , N z are the number of collocation points in the streamwise and spanwise di- 
rections, and a, f3 are the wavenumbers in these same two directions. The periods of the 
physical domain in the streamwise and spanwise directions are respectively 2tt / a and 27 r//?. 
In the absence of discontinuities in the solution, spectral collocation methods offer a high 
degree of accuracy and minimum dispersion and dissipation errors. Periodicity in the two 
horizontal directions allows the use of Fourier collocation methods. In the direction normal 
to the plate, Chebyshev collocation is the natural choice. The physical domain y £ [0,y maa; ] 
in the direction normal to the plate is mapped onto the computational domain 77 € [—1, l] 
through a convolution of two mappings. First, a hyperbolic tangent transformation takes 
77 to an intermediate variable ip according to 

, + Manh (i^) = ^. (10, 


In this equation, the hyperbolic tangent concentrates nodes about ipo, and A ip measures the 
width of the concentration region. t e measures the degree of influence of the hyperbolic 
tangent term. If t e = 0, there is no redistribution of nodes at this level. An algebraic 
stretching then maps ip onto the physical variable y £ [0 ,y max }. This mapping takes the 
form 


y _ yi/2ymaz(l + 

Umax lft{ymax 

One half the normal nodes are located between y = 0 and y = y\ji- With both mappings in 
hand, 770 , A 77 ipo and A ip are readily determined. The user specifies the following: l) t e , 2) 
3 / 0 ) the physical location about which the concentration occurs, 3) A y 0 , the concentration 
width in physical space, 4) y max , the maximum vertical extent of the physical domain, 
and 5) yi/ 2 - Admittedly, the selection is rather tedious, but it is only executed once to 
optimize the distribution of nodes relative to the initial perturbation wave profiles. The 
distribution of nodes is fixed during the simulation. Ay 0 and A ip are related via 


Ay 0 ( 12 ) 

V~Vo 

while ipo is computed from the algebraic mapping evaluated at y 0 . The two remaining 
unknowns, rjo and A 77 follow from the conditions 



vWU=-i = -1 
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= + 1 


(13) 


which leads to 


Arj = 


Vo = 


where 


T± = tanh 


2 + t e {T-+T+) 

<e(r + -r_) 

2 + * £ (r_ + t+) 

i ± 


Atjj 


(14) 


(15) 


All runs were performed with y max = 15, y x / 2 = 1, t ( = 0.8, y 0 = 1.2 and Ay 0 = 0.4. 


IV. Initial and Boundary Conditions 

In this paper, only the evolution of 2 dimensional waves are considered. Therefore, the 
initial conditions consist of a 2-D waves superimposed on a precalculated mean flow. Mean 
flow profiles are generated from the spectral solution to the similar compressible boundary- 
layer equations with zero pressure gradient and zero heat transfer at the wall [6], while the 
perturbation wave is a two-dimensional eigenfunction of the linearized compressible Navier- 
Stokes equations. The stability of flows is studied at a prescribed supercritical Reynolds 
number Re(x 0 ). After freezing the streamwise dependence of the mean flow profiles at 
x = xo, they are extended over the entire plate. Thus, the displacement thickness is 
constant, and the mean flow is a function only of the coordinate normal to the plate. 

The initial condition for u(x,y,z) is therefore 

u(x , y, z) = u m (y) + e 2D u 2D (y ) cos (0 2D (y) + ax) (16) 

with similar relations for the remaining 4 variables (v,w,p,p). The 2-D wave is an eigen- 
function to the linearized Navier-Stokes equations and is calculated with the spectral linear 
compressible stability code SPECLS [14]. The angle 0 2 d is the phase angle of the complex 
TS wave. The amplitudes u 2 d is normalized to unity. Therefore, e 2 r> is the amplitude of 
the two-dimensional perturbation wave with respect to the free-stream velocity. At high 
Mach numbers, e 2 r> cannot exceed a couple of percent lest the density and temperature 
maxima (which scale as M^) become excessive. 

Under the assumption of parallel flow, all the variables are periodic in the streamwise 
and spanwise directions. No-slip conditions are applied to the velocities at the wall which 
is adiabatic. In the far-field, the variables are either frozen at their initial values (which 
are vanishingly small) or calculated from zero stress conditions. These approaches both 
preclude the study of boundary-layer receptivity to incoming sound waves. The frozen 
boundary conditions also prevent waves from escaping into the free-stream. When the 
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far-field boundary is sufficiently removed from the wall, stability mechanisms which do not 
involve far-field interactions can still be successfully analyzed. 

Zero normal stress boundary conditions lead to an exponential decay of the solution 
in the free-stream, but not at the proper rate. However, it is plausible that the feedback 
from the free-stream into the interior of the domain is not as severe as for the frozen far- 
field conditions. These boundary conditions are derived from a variational formulation of 
the Navier-Stokes equations. The velocities and temperature are assumed to possess C 1 
continuity, while the pressure is only required to be C°. The natural boundary conditions 
of this formulation in the far-field are then [21] 


W x + U z = 0 

(17) 

Wy + Vg =0 

(18) 

- pRe + 2 w z + A(u z + v v + w z ) = 0 

(19) 

1 , dT, 

= -,Pr.R e M> dz ] 

(20) 


If ( u,v,w,p,T ) are deviations of ( u,v,w,p,T ) from free-stream values, the resulting lin- 
earized equations become 


u z — taw = 0 

(21) 

C2* 

N 

1 

C* 

e* 

ii 

o 

(22) 

— pRe + (2 + X'jWg + X(tau + — 0 

(23) 

W 7PrReMl Tz ~ ° 

(24) 

7 Ml,p z — p z + T Z = 0. 

(25) 


Derivatives in the vertical direction are evaluated using matrix multiplication. In this way, 
the unknowns (primitive variables at the last node) can be solved for simultaneously. The 
alternative are FFT’s, which lead to an iterative approach. The interior nodes have already 
been computed from the explicit time advancement. Therefore, the derivative of u at the 
far field (node N) is represented by 


du 

dy 


N 

= E D a u i 

V = Vi 


(26) 


where Z?, ; - is the collocation derivative matrix [2] and where the only unknown is u^. The 
system of equations (21)-(25) reduces to a linear system for the unknowns (un,vn, wn,Tn). 
Density has already been calculated through the continuity equation imposed at the outer 
boundary. A disadvantage comes from the higher truncation errors associated with matrix 
multiplication when compared to FFT’s for large number of collocation points. There 
are at least two remedies. First, an implicit, rather than an explicit time advancement 
scheme could incorporate these asymptotic boundary-conditions into the iterative inver- 
sion scheme. This could be done solely with FFT’s. The second approach is to adopt 
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a multidomain decomposition in the direction normal to the plate. An explicit scheme 
still requires matrix multiplication for derivative evaluation, but the number of nodes in 
each subdomain is small, and truncation errors will not affect the numerical results as 
strongly as in the single domain case. As the Mach number increases, the structure of 
the perturbation waves is mostly due to the density in a thin region centered about the 
critical layer whose distance away from the wall is approximately constant in terms of 8 * . 
Boundary conditions at the interface between two domains of a multidomain methods are 
continuity of all primitive variables (in the absence of shocks) and zero normal stress across 
the interface. A condition similar to zero normal stress can be derived from the energy 
equation. For best results, the interface conditions should be derived from the actual set 
of equations solved in the simulation. As yet, the multi-domain approach has not been 
implemented into the code. A more detailed description of this approach will be given 
elsewhere. All results presented below are based on the zero stress far-field conditions. 


V. Cost of Direct Simulation 

Timings on the Cray 2 at NASA Ames indicate that one iteration of the compressible 
Navier-Stokes equations is approximately twice the cost of one iteration of an incompress- 
ible spectral Navier-Stokes code. However, the ratio is actually higher for a given am- 
plification factor for the following reasons. Although the high Mach number regime does 
not require an implicit treatment of the acoustic terms, grid refinement in the vertical 
direction in the course of the simulation imposes a severe time step constraint due to the 
viscous terms. In the absence of grid stretching, a Chebyshev collocation method imposes 
a maximum explicit time step that scales as N ~ 4 , where N is the number of collocation 
points. With an algebraic grid stretching, points are more sparse in the free-stream, but 
concentrated at the origin, which decreases the time step further. Fortunately the hyper- 
bolic tangent stretching, which concentrates points near the critical layer at the expense 
of nodes near the wall, counters somewhat the effect of the algebraic mapping. The higher 
grid resolution required in the initial stages of the simulation due to high gradients in the 
density and temperature perturbation profiles in the critical layer region also increases the 
total CPU time. And last, but not least, compressibility effects weakens the instability. 
Second mode waves have a high spatial and temporal frequency. With the chosen parame- 
ters, the wave linearly grows by a factor of 1.06 every period. In the absence of non-linear 
effects, over 150 periods are required to achieve an amplification of e 9 , compared with less 
than 10 periods at low subsonic speeds. With an explicit code, simulation of 25 periods 
on a 8 X 2 x 65 grid requires approximately 10 CPU hours of Cray 2 at sustained speeds 
of 80 Mflops. At this stage, the growth of a single 2-D 2nd mode wave is still quasi-linear, 
and has only amplified by a factor 4. 
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A. Modal Analysis 


One of the many issues encountered in the treatment of compressible flows is the general- 
ization of modal energy, and intermodal energy transfer. These concepts are the backbone 
of incompressible transition theories [8]. Modal energy leads to the concept of energy 
cascade in turbulence [23], and permits the energy exchange processes between different 
modes to be computed. The standard definition of modal energy in incompressible flows 
is a consequence of the quadratic form of the kinetic energy. Let 

E = J v(x,y ) 2 dxdy = cst J \v{k,y )\ 2 dk dy (27) 

be the total kinetic energy in the system, where the second equality is a consequence of 
Parseval’s identity. v(x, y) and v(k, y) are Fourier Transforms of one another. The constant 
cst is inconsequential. The above definition of kinetic energy suggests that 

E k = cst J \v k (k,y)\ 2 dy (28) 

be interpreted as the integrated kinetic energy in mode k. Such a definition is no longer 
valid for compressible flows. The kinetic energy is now defined by the cubic nonlinearity 

E = \ J p(x,y)v(x,y ) 2 dxdy (29) 

which becomes in transformed space 

E — cst J p{k,y)v(l,y)v(—k — l,y) dkdl dy (30) 

which does not have a straightforward representation of the form 

E = J E(k) dk (31) 

One remedy, although not altogether justifiable, is to rewrite the kinetic energy as the 
product of velocity v and momentum m. For example, the kinetic energy of a photon must 
be defined this way because it is massless. If such a definition is accepted, the kinetic 
energy in mode k (intuitively) becomes 

E k = cst j (m(fc, y)v{-k, y) + c.c) dy (32) 

where the complex conjugate (c.c) is to insure that E k is real. It is then possible to 
determine contributions of any mode to any other mode. An equation describing the time 
evolution of E k can be derived, and intermodal transfer of kinetic energy determined. Of 
course it is by no means obvious that this is the physically appropriate definition for E, 
although it does reduce to the accepted form of the modal energy in the limit of zero Mach 
number. Another issue is: what fraction of internal energy is contained in each mode? The 
internal energy density is a quadratic function in p and T ; therefore, a definition similar 
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to that of Ek is possible. However, why should the internal energy be defined in terms of 
an energy density, while the kinetic energy is not? Ideally, the modal energy of different 
energy types should have consistent definitions. 

Experimentalists measure energy growth [10] (in the streamwise direction) in several 
ways, mostly with hot wire anemometry. The wire can be translated downstream at a 
constant height (y=const) above the plate. Alternatively, at each streamwise station, the 
y location can correspond to maximum growth of a specified variable. The growth rate 
can also be integrated in the y direction. Numerically, any one of these options is possible. 
We choose to calculate the instantaneous growth rate of the fundamental mode at each y 
location. Different variables have different growth rates and instantaneous growth rate of 
any variable can be measured. In this paper, only the growth of the streamwise velocity 
is considered because the presence of a critical layer might lead to interesting results. 


At any instant in time, say t*, u( k,y,t) has the representation 

U(k,y,t* +t) = u(k,y)e~ ia ^ v '^ t 
so that the growth rate is simply 


u 


*(v,o) = t- 


u 


t=0 


By comparing Eq. (34) with the growth rate computed from 

, i, u(k,y,t*+t) 

CT(y ’ i) = 7 log »(*,«, t-) ■ 


(33) 


(34) 


(35) 


the degree of linearity can be measured jus a function of t. If the growth of mode k is purely 
linear, both formulas will agree for all values of t. 


VI. Results 


In this paper, all results pertain to a fixed set of physical parameters. The free-stream 
Mach number is 4.5, = 110°I2, Re = 8000 (based on £*), a = 2.25 and (3 = 0.0. The 

choice of wavenumbers is motivated by the maximum temporal growth rate of 2-D second 
mode perturbation waves. The streamwise wavenumber produces the largest linear growth 
rate at the specified Reynolds number. 

The linear eigenfunctions were obtained from SPECLS [14] which solves the linearized 
Navier-Stokes equations with a spectral algorithm. Spectral initial conditions are neces- 
sary to prevent slight numerical oscillations in the mean flow and in the computed linear 
eigenfunction. Otherwise the non-linear calculation might be contaminated. Figures la-f 
compare the amplitude and phase plots of u, v, and T eigenfunctions of an unstable 2-D 
first mode wave (a = 0.6, f3 = 0, Re = 10000, Too = 110° 12 ) with the unstable second mode 
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defined above. The major differences between first and second mode disturbances are: 1) 
the first mode is viscous below a Mach number of about 3-4, and is nearly inviscid at 
higher Mach numbers, the second mode is always inviscid in nature; 2) the second mode 
phase velocity is close to unity; and 3) 3-D first modes are most unstable at angles between 
55° and 65° with respect to the mean flow, while 2-D second mode waves are most un- 
stable. Another feature worth noting for the particular choice of parameters, is the large 
vertical velocity component in the 2nd mode eigenfunction compared to that in the 1st 
mode wave. These differences might be explained by the factor 4 ratio of second to first 
mode streamwise wave numbers. As expected from Mack’s results, the second mode waves 
decay much faster than the 1st mode. The first mode waves decay on a viscous scale, while 
the 2nd mode decays on the inviscid scale. The 180° phase jumps in streamwise velocity 
occur in opposite directions (i.e. increasing towards the free-stream for the 2 nd mode, and 
decreasing for the 1 st mode). The significance of this, if any, is not known. 

The eigenfunctions calculated by the linear code are assumed to satisfy 

Upert [x, y, z ) = u' ( y) e *(«*- w 0 (36) 

so that the real part of u>, (cv r ) is the wave frequency and its imaginary part (u^) the growth 
rate. The wave is unstable if w,- is positive. The linear eigenvalue for the second mode 
perturbation used in the non-linear simulation is 

u> r = 2.046738 
cji = 2.149042 10" 2 

Based on the 2-D wave frequency, the temporal period is 3.07 in non-dimensional time 
units. From here on, time is always expressed in units in 2-D wave periods. Based on 
the above parameters, the streamwise phase velocity is c = 0.91, which corresponds to a 
critical layer located at y = 1.10. The mean flow profile also has a generalized inflection 
point * at y = 1.06. The mean temperature across the boundary layer varies by more than 
a factor of 4. This is in sharp contrast with the slowly varying mean temperature profiles 
at low subsonic Mach numbers. 

To test the code against linear theory, an initial eigenfunction with a 0.1% amplitude is 
superimposed on the mean flow and its evolution followed up to t = 3. The growth rate w,- 
and frequency (u; r ) are shown in figures 2a-b. The growth rate varies less than 0.4% over the 
interval 0 < y < 2. The frequency has a 0.1% variation in the same range of y. Integrated 
value of growth, when compared against prediction of linear theory, shows 4 decimal place 
agreement. In all the runs discussed, accuracy of the solution is measured from plots of the 
the Li norm of the Chebyshev series coefficients (integrated in the streamwise direction) 
versus the wavenumber. Grid resolution is chosen to insure that these spectra decrease by 
at least 8 orders of magnitude. 

Next, results from a run with a e 2D = 0.02 amplitude for streamwise perturbation 
*The generalized inflection point is the vertical coordinate y at which d/dy(p m (y)d/dy(U m (y))) = 0. 


10 



velocity, which is fairly large, are presented. The corresponding initial density fluctuation 
level is 9.2% 

One non-linear simulation of 25 periods in time was performed and the results are 
presented below. Initial streamwise velocity and density amplitude levels are 2% and 
9.2%. The remaining figures show the state of selected quantities at t=0,5,10,15 and 20 
(unit of time is one period). Figures 3a-b show the growth rate and frequency of the 
fundamental streamwise velocity, u 10 , labeled (1,0) for simplicity and to emphasize the 
two dimensionality of the problem. (In general, (k,0) refers to the kth mode, and (0,0) to 
the instantaneous mean flow.) The growth rate and the frequency are obtained from the 
coefficient of e l “* of the Fourier expansion of u calculated at two successive time levels. The 
major feature of interest is the sharp decrease in the local growth rate near the critical layer 
(y = 1.1). After a slight initial increase, the growth rate starts decreasing. If the growth 
rate reaches zero, the (1,0) mode will have reached a saturated state. The frequency, and 
therefore the phase velocity is decreasing slightly towards zero, albeit at a slower rate than 
the growth rate (fig. 3b) . The growth rate and frequency show a marked change near the 
wall at period 0. There are several possible reasons for this, but the phenomena has not 
been looked into in detail. Typically, numerically generated growth rates are obtained from 
the integrated energy, which wipes out any oscillations due to inconsistent or inaccurate 
initial conditions. The mean flow and the initial profiles were obtained spectrally, with two 
different codes and grids, and spectral interpolation was performed to transfer information 
from one grid to another. Because the mean flow was computed on a smaller grid than 
were the fundamental profiles, the interpolation generated slight oscillations in the second 
derivative of the mean flow. This in turn affects the local growth profiles. It is obvious 
however from the figure that the initial slight inconsistencies are washed away at later 
times. Note the sharp peak near the critical layer. To understand better the origin of 
this critical layer structure, Fourier transform the streamwse momentum equation in the x 
direction, and only consider the terms that contribute to the fundamental mode. Neglecting 
the viscous terms, this leads to the equation 

p 0 u u + ta(/o 0 uoUi + pi) + p 0 (viUo v + v 0 ui v )Lap-iu\ + p_iViUi„ + Pi{viu. ly + v_iUi v ). (37) 

if u\ is assumed to vary as e -UJt (where u is a function of y and t), the instantaneous 
growth rate and frequency are obtained through after dividing this equation by p 0 «i: 

— loj = linear terms + cubic terms (38) 

where 

linear terms = ta(u 0 + — ) + — ■ llQy — v o u iy (39) 

Ul U! ’ 

and 

cubic terms = ± (40) 

Po PoUl 

The growth rate is the real part of Eq. (38). The linear terms contribute a zeroth order 
growth rate, while the growth rate due to the cubic nonlinearities is quadratic in the 
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amplitude of the fundamental mode. In figures 4a-b, the growth rate u r is plotted as a 
function of the normal coordinate y. Figure 4a plots the combination of linear terms, 
excluding the viscous interaction, which explains the sharp increase in growth near the 
origin and the structure at the critical layer, which are removed by the viscous effects. If 
the viscous effects were included, and the wave amplitude sufficiently small, figs. 4a and 
2a would be identical. Cubic nonlinearities are plotted in figure 4b. One notes the absence 
of growth rate everywhere, except in the critical layer and in a thin layer near the wall. 
The nonlinearities near the wall can be explained by consideration of the viscous terms. 
The viscosity is a function of temperature and becomes, after taylor expansion, 

n{T) » Mo + MorTi + \T-iTiHqtt + fTLiTjVorrr- (41) 

Derivatives of the mean viscosity p 0 with respect to the mean temperature are denoted 
by the subscripts T. Since the viscous stress terms are quadratic in the viscosity and 
the primitive variables, only the mean and fundamental components of viscosity need be 
retained. The terms that contribute to the linear equations (in the fundamental amplitude) 
are 

Mo + MotTx (42) 

while 

\T_xTiHott + \T-iTIhqttt (43) 

contributes to the cubic nonlinearities. The cubic nonlinearities drop out if p, is independent 
of temperature. More detailed mapping of these various terms is under progress. Figures 
5-7 show the fundamental, first harmonic, and distorted mean profiles. A comparison of fig. 
5a and fig. 5c shows that Ui and pi grow at different rates. Initially, ui has a normalized 
maximum of 0.01, and p has its maximum at 0.046, i.e. a pi to ui ratio of 4.6. At t = 20, 
this ratio is 2.62, indicating that the peak u\ grows faster than the peak p\. Of course 
these growth rates may be governed by different processes. Near the wall, viscous effects 
dominate, while p has its maximum at the critical layer, where the cubic interactions 
dominate. The mean vertical flow is plotted in fig. 6 and has two strong peaks. Near 
the wall, the peak is broad, compared to the narrow peak near the critical layer. These 
structures are best explained by a perturbation methods approach, which is currently under 
investigation. Figures 7a^c show the first harmonic components of u,v,p. The maximum 
amplitude of these waves is approximately an order of magnitude smaller than that of the 
fundamental. This suggests that a weakly non-linear asymptotic theory is sufficient to 
predict the evolution of a single two-dimensional mode, perhaps including the structure of 
its saturated state. To support this hypothesis, a second calculation was performed with 
the same initial conditions. However, all streamwise modes above the fundamental were set 
to zero after each iteration s . Growth rate curves at t = 10 are virtually identical to that 
shown in fig. 3. This indicates that the higher streamwise modes only play a secondary 
role in the evolution towards a saturated state. As shown above, the leading role is that 
due to the cubic nonlinearities in the momentum equation A perturbation analysis should 
allow a large number of parametric studies relatively cheaply. Contributions to the growth 

More precisely they were set to zero at each stage of the 3rd order Runga-Kutta algorithm 
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rate from various sources will be determined, and the simulation carried out further than 
otherwise possible. 


Finally, it is interesting to understand the mechanisms that pertain to vorticity pro- 
duction. The compressible analogue to the incompressible vorticity equation is 

^ + v • Vw = 3 • Vv — toV-v — Vp -1 x Vp + 2 pVx 


dt 


P _1 V-(r-|V.t7/) (44) 


where 3 = V x v. Terms not present in the limit of incompressibility are the dilatation term 
wV-v and the baroclinic torque Vp -1 x Vp. The vorticity production due to dilatation 
and baroclinic torque (average over streamwise direction) are shown in figs 8a-b. Figure 8c 
shows the combined effects of dilation, convection and baroclinic torque on the production 
of vorticity. In all 3 plots, the production of vorticity is localized in the near wall region, 
and at the critical layer. The rms perturbation vorticity is shown in fig. 9 after 20 periods 
(after subtraction of the mean flow vorticity) . 


VII. Conclusions 


A three-dimensional, fully spectral code developed for the direct simulation of the wall- 
bounded compressible flows has been applied to the nonlinear evolution of a single unstable 
second mode wave. After an slight initial growth, the growth rate of this fundamental wave 
decreases in time. This suggest that the second mode is evolving towards a saturated state. 
Perturbation methods are sometimes well-suited for the study of such states. To ascertain 
whether such an approach is called for, all streamwise modes above the fundamental were 
continuously set to zero at every time step. This led to a time evolution of the fundamental 
waves almost identical with the case where all modes were allowed to interact with each 
other. This agreement was also found for the structure of the time-dependent instanta- 
neous growth rate of u, even in the critical layer region. The conclusion is that the cubic 
interaction terms that only involve the fundamental waves are probably responsible for the 
nonlinearities in the critical layer, and that the density fluctuations play a fundamental 
role. 

A final note. The initial perturbation wave was chosen to maximize its instability. 
Nayfeh [20] found that stable 2-D second mode waves were susceptible to strong secondary 
instabilities. It is therefore still not clear which are the actual waves that will be most 
responsible for the final transition towards turbulence. These issues will be addressed in 
the future. 
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Figure lc 
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Figure Id 
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Figure 2b 
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Figure 3a 
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Figure 3b 
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Figure 4a 
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Figure 4b 



2nd mode run19 M=4.5 



o 


y- 

y- 



T- 

1— 


o 

o 

o 

o 

o 

o 

o 

1 

I 

1 

1 

1 

1 

1 

0 

0 

* 

• 

• 

0 

• 

to 

o 

10 

o 

10 

o 

o 

n 

ro 

CM 

M 

r- 


in 

6 

6 

6 

6 

d 

6 

o 


(OL) |duuD n 


28 


Figure 5a 
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Figure 5c 
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Figure 7a 
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Figure 7b 
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Figure 7c 







2nd mod© run19 M=4.5 



35 


Figure 8a 
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Figure 8b 
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